library(rstan)
## Loading required package: ggplot2
## Loading required package: StanHeaders
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
library(ggmcmc)
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: tidyr
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:rstan':
## 
##     extract
library(bayesplot)
## This is bayesplot version 1.6.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
library(bridgesampling)
load("data_for_stan.RData")
list_for_stan <- list(nY = nrow(data_for_stan), nS = max(data_for_stan$participant), Subj = data_for_stan$participant, size_diff = data_for_stan$size, num_of_trials = data_for_stan$number_of_fixational_trials, Y = data_for_stan$answer.keys)

Оценка связи типа иллюзии и разницы между стимулами

fit_1 <- stan("bayes_model_1.stan", data = list_for_stan, iter = 6000, warmup = 2000)
## 
## SAMPLING FOR MODEL 'bayes_model_1' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000102 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.02 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 1: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 1: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 1: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 1: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 1: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 1: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 1: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 1: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.01642 seconds (Warm-up)
## Chain 1:                2.49926 seconds (Sampling)
## Chain 1:                3.51568 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'bayes_model_1' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 3.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 2: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 2: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 2: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 2: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 2: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 2: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 2: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 2: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.965543 seconds (Warm-up)
## Chain 2:                2.53149 seconds (Sampling)
## Chain 2:                3.49704 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'bayes_model_1' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 3.8e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 3: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 3: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 3: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 3: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 3: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 3: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 3: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 3: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.03701 seconds (Warm-up)
## Chain 3:                2.29807 seconds (Sampling)
## Chain 3:                3.33508 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'bayes_model_1' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 4.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 4: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 4: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 4: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 4: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 4: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 4: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 4: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 4: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.06415 seconds (Warm-up)
## Chain 4:                2.57186 seconds (Sampling)
## Chain 4:                3.636 seconds (Total)
## Chain 4:
print(fit_1)
## Inference for Stan model: bayes_model_1.
## 4 chains, each with iter=6000; warmup=2000; thin=1; 
## post-warmup draws per chain=4000, total post-warmup draws=16000.
## 
##              mean se_mean   sd    2.5%     25%     50%     75%   97.5%
## intercept    1.06    0.01 0.43    0.24    0.77    1.05    1.35    1.93
## beta_mu     -0.94    0.00 0.32   -1.60   -1.14   -0.92   -0.72   -0.34
## beta_sd      1.12    0.00 0.27    0.70    0.93    1.09    1.28    1.76
## beta[1]     -0.48    0.00 0.34   -1.15   -0.71   -0.48   -0.25    0.19
## beta[2]     -0.06    0.00 0.48   -0.97   -0.39   -0.08    0.24    0.95
## beta[3]     -1.75    0.01 0.68   -3.32   -2.14   -1.67   -1.28   -0.67
## beta[4]     -2.40    0.01 0.85   -4.34   -2.89   -2.29   -1.78   -1.06
## beta[5]     -1.33    0.01 0.52   -2.47   -1.65   -1.28   -0.96   -0.42
## beta[6]      0.39    0.01 0.58   -0.56    0.00    0.34    0.72    1.70
## beta[7]     -0.01    0.00 0.38   -0.74   -0.27   -0.03    0.23    0.76
## beta[8]     -0.04    0.00 0.36   -0.70   -0.28   -0.05    0.20    0.70
## beta[9]     -0.38    0.00 0.40   -1.16   -0.65   -0.38   -0.11    0.43
## beta[10]    -2.13    0.01 0.85   -4.09   -2.60   -2.01   -1.53   -0.81
## beta[11]    -1.49    0.01 0.46   -2.49   -1.77   -1.45   -1.18   -0.69
## beta[12]    -1.46    0.01 0.50   -2.52   -1.76   -1.42   -1.11   -0.59
## beta[13]    -1.69    0.01 0.90   -3.71   -2.19   -1.59   -1.06   -0.26
## beta[14]    -1.09    0.00 0.40   -1.94   -1.34   -1.06   -0.81   -0.38
## beta[15]    -2.33    0.01 0.83   -4.27   -2.79   -2.22   -1.75   -1.05
## beta[16]    -0.88    0.00 0.36   -1.62   -1.10   -0.86   -0.63   -0.22
## beta[17]     0.09    0.00 0.36   -0.58   -0.16    0.07    0.32    0.84
## beta[18]    -1.33    0.01 0.74   -2.99   -1.77   -1.26   -0.81   -0.06
## beta[19]    -2.46    0.01 0.82   -4.37   -2.91   -2.35   -1.89   -1.20
## beta[20]    -0.16    0.00 0.34   -0.82   -0.39   -0.17    0.06    0.54
## beta[21]    -1.63    0.01 0.88   -3.68   -2.12   -1.52   -1.02   -0.23
## beta[22]    -1.30    0.00 0.43   -2.26   -1.57   -1.27   -1.00   -0.54
## beta[23]     0.25    0.00 0.35   -0.38    0.00    0.23    0.47    1.00
## beta[24]     0.58    0.00 0.47   -0.22    0.25    0.54    0.86    1.64
## beta[25]    -0.37    0.00 0.47   -1.30   -0.68   -0.37   -0.06    0.55
## beta[26]    -2.15    0.01 0.87   -4.14   -2.63   -2.04   -1.54   -0.76
## lp__      -124.29    0.08 4.85 -134.92 -127.28 -123.90 -120.90 -115.92
##           n_eff Rhat
## intercept  4215    1
## beta_mu    5378    1
## beta_sd    4755    1
## beta[1]    8271    1
## beta[2]   12277    1
## beta[3]    8287    1
## beta[4]    7443    1
## beta[5]    9931    1
## beta[6]   11904    1
## beta[7]    9502    1
## beta[8]    9462    1
## beta[9]   11787    1
## beta[10]   7840    1
## beta[11]   8144    1
## beta[12]   8523    1
## beta[13]  10069    1
## beta[14]   9490    1
## beta[15]   8002    1
## beta[16]   9101    1
## beta[17]   9959    1
## beta[18]  11723    1
## beta[19]   7063    1
## beta[20]  10318    1
## beta[21]  10001    1
## beta[22]   8502    1
## beta[23]   9879    1
## beta[24]  10404    1
## beta[25]  12476    1
## beta[26]   7888    1
## lp__       3422    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Jan  4 16:42:49 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
posterior_1 <- as.array(fit_1)
mcmc_areas(posterior_1, regex_pars = c("beta"), prob = 0.8, prob_outer = 0.95, point_est = "mean")

Диагностика

posterior_1b <- ggs(fit_1)
ggmcmc(D = posterior_1b, file = NULL, plot = 'ggs_traceplot') 
## Plotting traceplots

## Time taken to generate the report: 16 seconds.
ggmcmc(D = posterior_1b, file = NULL, plot = 'ggs_compare_partial') 
## Plotting comparison of partial and full chain

## Time taken to generate the report: 19 seconds.
ggmcmc(D = posterior_1b, file = NULL, plot = 'ggs_autocorrelation') 
## Plotting autocorrelation plots

## Time taken to generate the report: 15 seconds.

Модель с индивидуальными смещениями

fit_2 <- stan("bayes_model_2.stan", data = list_for_stan, iter = 6000, warmup = 2000)
## 
## SAMPLING FOR MODEL 'bayes_model_2' NOW (CHAIN 1).
## Chain 1: Rejecting initial value:
## Chain 1:   Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: Random variable is -1.45713, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000139 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.39 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 1: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 1: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 1: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 1: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 1: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 1: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 1: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 1: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.04328 seconds (Warm-up)
## Chain 1:                1.8008 seconds (Sampling)
## Chain 1:                2.84408 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'bayes_model_2' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 2: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 2: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 2: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 2: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 2: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 2: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 2: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 2: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.15235 seconds (Warm-up)
## Chain 2:                1.67036 seconds (Sampling)
## Chain 2:                2.82271 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'bayes_model_2' NOW (CHAIN 3).
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: beta_lpdf: Random variable is 1.19705, but must be less than or equal to 1  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: beta_lpdf: Random variable is -0.336032, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: beta_lpdf: Random variable is 1.52124, but must be less than or equal to 1  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: beta_lpdf: Random variable is -0.0588813, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: beta_lpdf: Random variable is -0.963093, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: beta_lpdf: Random variable is 1.5279, but must be less than or equal to 1  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 3: 
## Chain 3: Gradient evaluation took 4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.4 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 3: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 3: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 3: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 3: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 3: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 3: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 3: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 3: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.26678 seconds (Warm-up)
## Chain 3:                2.88996 seconds (Sampling)
## Chain 3:                4.15674 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'bayes_model_2' NOW (CHAIN 4).
## Chain 4: Rejecting initial value:
## Chain 4:   Error evaluating the log probability at the initial value.
## Chain 4: Exception: beta_lpdf: Random variable is -0.264047, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 4: Rejecting initial value:
## Chain 4:   Error evaluating the log probability at the initial value.
## Chain 4: Exception: beta_lpdf: Random variable is 1.2442, but must be less than or equal to 1  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 4: Rejecting initial value:
## Chain 4:   Error evaluating the log probability at the initial value.
## Chain 4: Exception: beta_lpdf: Random variable is -1.0056, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 4: Rejecting initial value:
## Chain 4:   Error evaluating the log probability at the initial value.
## Chain 4: Exception: beta_lpdf: Random variable is -0.0367961, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 4: Rejecting initial value:
## Chain 4:   Error evaluating the log probability at the initial value.
## Chain 4: Exception: beta_lpdf: Random variable is -1.9334, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 4: Rejecting initial value:
## Chain 4:   Error evaluating the log probability at the initial value.
## Chain 4: Exception: beta_lpdf: Random variable is -1.09291, but must be >= 0!  (in 'model47ca61cf6d02_bayes_model_2' at line 20)
## 
## Chain 4: 
## Chain 4: Gradient evaluation took 6.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 4: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 4: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 4: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 4: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 4: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 4: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 4: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 4: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.21527 seconds (Warm-up)
## Chain 4:                1.7629 seconds (Sampling)
## Chain 4:                2.97817 seconds (Total)
## Chain 4:
## Warning: There were 9082 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
print(fit_2)
## Inference for Stan model: bayes_model_2.
## 4 chains, each with iter=6000; warmup=2000; thin=1; 
## post-warmup draws per chain=4000, total post-warmup draws=16000.
## 
##                  mean se_mean   sd    2.5%     25%     50%     75%   97.5%
## intercept_mu     0.47    0.01 0.28    0.03    0.24    0.46    0.70    0.97
## intercept_sd     2.05    0.01 0.46    1.32    1.72    1.99    2.31    3.10
## intercept[1]     1.52    0.01 0.65    0.26    1.08    1.52    1.96    2.81
## intercept[2]     2.54    0.02 0.99    0.79    1.86    2.48    3.15    4.67
## intercept[3]    -0.67    0.01 0.98   -2.72   -1.29   -0.61    0.02    1.07
## intercept[4]    -2.20    0.02 1.34   -5.23   -2.95   -2.02   -1.26   -0.07
## intercept[5]    -0.18    0.01 0.82   -1.91   -0.70   -0.15    0.36    1.38
## intercept[6]     2.48    0.02 1.02    0.64    1.79    2.41    3.11    4.64
## intercept[7]     2.23    0.01 0.71    0.92    1.74    2.21    2.68    3.75
## intercept[8]     2.11    0.01 0.76    0.67    1.58    2.09    2.61    3.64
## intercept[9]     1.94    0.01 0.87    0.32    1.35    1.90    2.50    3.72
## intercept[10]   -1.79    0.02 1.41   -5.04   -2.57   -1.64   -0.81    0.52
## intercept[11]   -0.82    0.01 0.78   -2.48   -1.31   -0.78   -0.29    0.59
## intercept[12]    0.06    0.01 0.65   -1.26   -0.37    0.08    0.51    1.31
## intercept[13]   -0.91    0.02 1.59   -4.46   -1.88   -0.78    0.20    1.84
## intercept[14]   -0.29    0.01 0.83   -2.07   -0.81   -0.24    0.26    1.25
## intercept[15]   -2.15    0.02 1.33   -5.23   -2.88   -1.98   -1.22    0.01
## intercept[16]    0.23    0.01 0.75   -1.29   -0.26    0.25    0.75    1.66
## intercept[17]    2.62    0.02 0.83    1.09    2.05    2.59    3.14    4.34
## intercept[18]    0.32    0.02 1.18   -2.10   -0.44    0.35    1.12    2.53
## intercept[19]   -2.42    0.02 1.30   -5.52   -3.12   -2.24   -1.50   -0.34
## intercept[20]    1.61    0.01 0.84   -0.01    1.04    1.61    2.19    3.28
## intercept[21]   -0.76    0.03 1.68   -4.59   -1.72   -0.62    0.39    2.11
## intercept[22]    0.09    0.01 0.65   -1.19   -0.34    0.11    0.53    1.34
## intercept[23]    3.01    0.02 0.82    1.52    2.45    2.96    3.55    4.77
## intercept[24]    2.90    0.01 0.81    1.41    2.34    2.86    3.41    4.61
## intercept[25]    1.41    0.01 0.94   -0.40    0.79    1.39    2.02    3.31
## intercept[26]   -1.78    0.03 1.47   -5.12   -2.57   -1.60   -0.78    0.55
## beta            -0.57    0.00 0.17   -0.90   -0.69   -0.58   -0.46   -0.24
## lp__          -142.22    0.09 4.50 -151.99 -145.01 -141.90 -139.01 -134.43
##               n_eff Rhat
## intercept_mu   2148    1
## intercept_sd   2793    1
## intercept[1]   3557    1
## intercept[2]   3477    1
## intercept[3]   4954    1
## intercept[4]   4627    1
## intercept[5]   4829    1
## intercept[6]   4304    1
## intercept[7]   3547    1
## intercept[8]   2925    1
## intercept[9]   3974    1
## intercept[10]  4124    1
## intercept[11]  3553    1
## intercept[12]  3857    1
## intercept[13]  5360    1
## intercept[14]  3673    1
## intercept[15]  4555    1
## intercept[16]  3396    1
## intercept[17]  3053    1
## intercept[18]  6034    1
## intercept[19]  4446    1
## intercept[20]  3712    1
## intercept[21]  4482    1
## intercept[22]  3585    1
## intercept[23]  2947    1
## intercept[24]  3091    1
## intercept[25]  4606    1
## intercept[26]  3256    1
## beta           1500    1
## lp__           2751    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Jan  4 16:47:12 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
posterior_2 <- as.array(fit_2)
mcmc_areas(posterior_2, regex_pars = c("intercept"), prob = 0.8, prob_outer = 0.95, point_est = "mean")

Диагностика

posterior_2b <- ggs(fit_2)
ggmcmc(D = posterior_2b, file = NULL, plot = 'ggs_traceplot') 
## Plotting traceplots

## Time taken to generate the report: 14 seconds.
ggmcmc(D = posterior_2b, file = NULL, plot = 'ggs_compare_partial') 
## Plotting comparison of partial and full chain

## Time taken to generate the report: 18 seconds.
ggmcmc(D = posterior_2b, file = NULL, plot = 'ggs_autocorrelation') 
## Plotting autocorrelation plots

## Time taken to generate the report: 15 seconds.

Модель с индивидуальными смещениями и индивидуальными размерами эффекта

fit_3 <- stan("bayes_model_3.stan", data = list_for_stan, iter = 6000, warmup = 2000)
## 
## SAMPLING FOR MODEL 'bayes_model_3' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.56 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 1: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 1: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 1: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 1: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 1: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 1: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 1: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 1: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.91152 seconds (Warm-up)
## Chain 1:                2.86213 seconds (Sampling)
## Chain 1:                4.77364 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'bayes_model_3' NOW (CHAIN 2).
## Chain 2: Rejecting initial value:
## Chain 2:   Error evaluating the log probability at the initial value.
## Chain 2: Exception: beta_lpdf: Random variable is 1.16891, but must be less than or equal to 1  (in 'model47ca2e7c7cdd_bayes_model_3' at line 21)
## 
## Chain 2: Rejecting initial value:
## Chain 2:   Error evaluating the log probability at the initial value.
## Chain 2: Exception: beta_lpdf: Random variable is 1.61189, but must be less than or equal to 1  (in 'model47ca2e7c7cdd_bayes_model_3' at line 21)
## 
## Chain 2: Rejecting initial value:
## Chain 2:   Error evaluating the log probability at the initial value.
## Chain 2: Exception: beta_lpdf: Random variable is -0.947509, but must be >= 0!  (in 'model47ca2e7c7cdd_bayes_model_3' at line 21)
## 
## Chain 2: Rejecting initial value:
## Chain 2:   Error evaluating the log probability at the initial value.
## Chain 2: Exception: beta_lpdf: Random variable is -1.26298, but must be >= 0!  (in 'model47ca2e7c7cdd_bayes_model_3' at line 21)
## 
## Chain 2: 
## Chain 2: Gradient evaluation took 4.6e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.46 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 2: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 2: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 2: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 2: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 2: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 2: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 2: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 2: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.86638 seconds (Warm-up)
## Chain 2:                3.48361 seconds (Sampling)
## Chain 2:                5.34999 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'bayes_model_3' NOW (CHAIN 3).
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: beta_lpdf: Random variable is 1.24996, but must be less than or equal to 1  (in 'model47ca2e7c7cdd_bayes_model_3' at line 21)
## 
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: beta_lpdf: Random variable is 1.87894, but must be less than or equal to 1  (in 'model47ca2e7c7cdd_bayes_model_3' at line 21)
## 
## Chain 3: 
## Chain 3: Gradient evaluation took 4.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 3: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 3: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 3: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 3: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 3: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 3: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 3: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 3: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.20559 seconds (Warm-up)
## Chain 3:                3.05649 seconds (Sampling)
## Chain 3:                5.26208 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'bayes_model_3' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 4.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 6000 [  0%]  (Warmup)
## Chain 4: Iteration:  600 / 6000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1200 / 6000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1800 / 6000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2001 / 6000 [ 33%]  (Sampling)
## Chain 4: Iteration: 2600 / 6000 [ 43%]  (Sampling)
## Chain 4: Iteration: 3200 / 6000 [ 53%]  (Sampling)
## Chain 4: Iteration: 3800 / 6000 [ 63%]  (Sampling)
## Chain 4: Iteration: 4400 / 6000 [ 73%]  (Sampling)
## Chain 4: Iteration: 5000 / 6000 [ 83%]  (Sampling)
## Chain 4: Iteration: 5600 / 6000 [ 93%]  (Sampling)
## Chain 4: Iteration: 6000 / 6000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.05636 seconds (Warm-up)
## Chain 4:                20.156 seconds (Sampling)
## Chain 4:                22.2123 seconds (Total)
## Chain 4:
## Warning: There were 10422 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
print(fit_3)
## Inference for Stan model: bayes_model_3.
## 4 chains, each with iter=6000; warmup=2000; thin=1; 
## post-warmup draws per chain=4000, total post-warmup draws=16000.
## 
##                  mean se_mean    sd    2.5%     25%     50%     75%
## intercept_mu     0.58    0.01  0.27    0.05    0.37    0.61    0.82
## intercept_sd     1.35    0.03  0.62    0.25    0.89    1.32    1.75
## intercept[1]     1.65    0.02  1.02   -0.06    0.92    1.55    2.28
## intercept[2]     1.92    0.03  1.23    0.00    1.03    1.76    2.67
## intercept[3]     0.40    0.02  1.10   -1.89   -0.23    0.46    1.07
## intercept[4]    -0.62    0.04  1.40   -3.73   -1.44   -0.42    0.39
## intercept[5]     0.45    0.02  0.99   -1.61   -0.15    0.49    1.07
## intercept[6]     1.08    0.02  1.12   -1.00    0.35    0.98    1.75
## intercept[7]     1.59    0.03  1.01   -0.11    0.89    1.49    2.21
## intercept[8]     1.08    0.02  1.01   -0.77    0.43    0.98    1.66
## intercept[9]     1.79    0.03  1.19   -0.12    0.91    1.65    2.51
## intercept[10]   -0.35    0.03  1.39   -3.63   -1.12   -0.15    0.60
## intercept[11]   -0.12    0.02  1.04   -2.42   -0.75    0.00    0.60
## intercept[12]    1.14    0.02  0.95   -0.55    0.51    1.05    1.70
## intercept[13]    0.07    0.03  1.43   -3.25   -0.66    0.27    0.95
## intercept[14]   -0.06    0.02  1.09   -2.50   -0.70    0.06    0.67
## intercept[15]   -0.58    0.04  1.45   -4.03   -1.39   -0.35    0.46
## intercept[16]    0.21    0.02  1.04   -1.99   -0.41    0.29    0.87
## intercept[17]    1.61    0.03  1.15   -0.32    0.80    1.48    2.30
## intercept[18]    0.82    0.02  1.20   -1.62    0.13    0.78    1.50
## intercept[19]   -0.64    0.04  1.42   -3.94   -1.46   -0.42    0.41
## intercept[20]    0.30    0.02  1.17   -2.22   -0.36    0.40    1.01
## intercept[21]    0.17    0.02  1.44   -3.05   -0.58    0.32    1.02
## intercept[22]    1.01    0.02  0.92   -0.73    0.42    0.96    1.56
## intercept[23]    1.73    0.03  1.15   -0.17    0.90    1.60    2.45
## intercept[24]    0.76    0.02  1.09   -1.41    0.10    0.73    1.41
## intercept[25]    0.96    0.02  1.12   -1.18    0.26    0.89    1.62
## intercept[26]   -0.38    0.03  1.37   -3.64   -1.13   -0.18    0.56
## beta_mu         -0.74    0.01  0.27   -1.31   -0.90   -0.73   -0.56
## beta_sd          0.93    0.01  0.32    0.38    0.72    0.91    1.12
## beta[1]         -0.71    0.01  0.50   -1.76   -1.01   -0.68   -0.37
## beta[2]         -0.37    0.01  0.59   -1.57   -0.74   -0.36    0.01
## beta[3]         -1.38    0.02  0.71   -2.97   -1.78   -1.31   -0.89
## beta[4]         -1.63    0.02  0.90   -3.65   -2.14   -1.54   -1.00
## beta[5]         -1.02    0.01  0.58   -2.27   -1.37   -0.99   -0.63
## beta[6]          0.33    0.01  0.66   -0.81   -0.13    0.28    0.71
## beta[7]         -0.23    0.01  0.53   -1.29   -0.58   -0.24    0.10
## beta[8]         -0.05    0.01  0.49   -1.01   -0.37   -0.05    0.27
## beta[9]         -0.61    0.01  0.54   -1.76   -0.94   -0.59   -0.24
## beta[10]        -1.51    0.02  0.88   -3.59   -1.99   -1.40   -0.92
## beta[11]        -0.95    0.01  0.55   -2.05   -1.31   -0.94   -0.59
## beta[12]        -1.42    0.02  0.65   -2.88   -1.80   -1.36   -0.97
## beta[13]        -1.26    0.02  0.84   -3.24   -1.73   -1.15   -0.69
## beta[14]        -0.65    0.01  0.51   -1.66   -0.98   -0.64   -0.32
## beta[15]        -1.60    0.02  0.85   -3.52   -2.09   -1.52   -1.01
## beta[16]        -0.55    0.01  0.49   -1.52   -0.87   -0.55   -0.24
## beta[17]        -0.12    0.01  0.51   -1.16   -0.45   -0.12    0.21
## beta[18]        -1.10    0.01  0.73   -2.75   -1.53   -1.05   -0.61
## beta[19]        -1.69    0.02  0.86   -3.61   -2.18   -1.61   -1.09
## beta[20]         0.07    0.01  0.49   -0.84   -0.26    0.05    0.37
## beta[21]        -1.24    0.02  0.80   -3.11   -1.69   -1.13   -0.71
## beta[22]        -1.22    0.01  0.54   -2.42   -1.55   -1.19   -0.85
## beta[23]        -0.01    0.01  0.51   -1.01   -0.35    0.00    0.32
## beta[24]         0.65    0.01  0.64   -0.52    0.22    0.60    1.02
## beta[25]        -0.34    0.01  0.59   -1.52   -0.72   -0.34    0.05
## beta[26]        -1.49    0.02  0.87   -3.47   -1.96   -1.39   -0.89
## lp__          -132.11    0.71 13.23 -155.09 -140.49 -133.27 -125.18
##                 97.5% n_eff Rhat
## intercept_mu     0.98  1552    1
## intercept_sd     2.69   541    1
## intercept[1]     3.85  1729    1
## intercept[2]     4.73  1362    1
## intercept[3]     2.52  3066    1
## intercept[4]     1.60  1540    1
## intercept[5]     2.40  2635    1
## intercept[6]     3.53  2431    1
## intercept[7]     3.82  1616    1
## intercept[8]     3.27  2958    1
## intercept[9]     4.53  1355    1
## intercept[10]    1.86  1983    1
## intercept[11]    1.65  1800    1
## intercept[12]    3.20  1986    1
## intercept[13]    2.54  2303    1
## intercept[14]    1.86  2201    1
## intercept[15]    1.59  1478    1
## intercept[16]    2.18  2575    1
## intercept[17]    4.19  1643    1
## intercept[18]    3.39  2712    1
## intercept[19]    1.56  1464    1
## intercept[20]    2.51  2690    1
## intercept[21]    2.78  3428    1
## intercept[22]    2.96  2273    1
## intercept[23]    4.23  1554    1
## intercept[24]    3.07  2371    1
## intercept[25]    3.32  3240    1
## intercept[26]    1.83  1559    1
## beta_mu         -0.26  1408    1
## beta_sd          1.64  1252    1
## beta[1]          0.20  2092    1
## beta[2]          0.81  2150    1
## beta[3]         -0.20  2208    1
## beta[4]         -0.15  1925    1
## beta[5]          0.05  2739    1
## beta[6]          1.78  2232    1
## beta[7]          0.80  1536    1
## beta[8]          0.91  3057    1
## beta[9]          0.39  1718    1
## beta[10]        -0.12  2157    1
## beta[11]         0.11  2255    1
## beta[12]        -0.28  1467    1
## beta[13]         0.13  2355    1
## beta[14]         0.38  2498    1
## beta[15]        -0.19  2209    1
## beta[16]         0.42  2352    1
## beta[17]         0.87  2153    1
## beta[18]         0.18  2840    1
## beta[19]        -0.25  1630    1
## beta[20]         1.10  2810    1
## beta[21]         0.03  2739    1
## beta[22]        -0.29  1942    1
## beta[23]         1.00  1950    1
## beta[24]         2.05  1988    1
## beta[25]         0.84  3469    1
## beta[26]        -0.08  2045    1
## lp__          -101.60   346    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Jan  4 16:50:44 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
posterior_3 <- as.array(fit_3)
mcmc_areas(posterior_3, regex_pars = c("beta"), prob = 0.8, prob_outer = 0.95, point_est = "mean")

Диагностика

posterior_3b <- ggs(fit_3)
ggmcmc(D = posterior_3b, file = NULL, plot = 'ggs_traceplot') 
## Plotting traceplots

## Time taken to generate the report: 26 seconds.
ggmcmc(D = posterior_3b, file = NULL, plot = 'ggs_compare_partial') 
## Plotting comparison of partial and full chain

## Time taken to generate the report: 40 seconds.
ggmcmc(D = posterior_3b, file = NULL, plot = 'ggs_autocorrelation') 
## Plotting autocorrelation plots

## Time taken to generate the report: 31 seconds.